/* 
 * File:   multisurfalign_grid3d.h
 * Author: boyan
 *
 * Created on April 17, 2009, 12:59 PM
 */

#ifndef _MULTISURFALIGN_GRID3D_H
#define	_MULTISURFALIGN_GRID3D_H

#include "SurfAlign_grid3d.h"

#include <list>
#include <vector>
#include <set>
#include <bitset>
#include <iomanip>
#include <algorithm>

class MultiSurfAlign_grid3d : public SurfAlign_grid3d
{
public:

    MultiSurfAlign_grid3d(grid3d_t *g);
    MultiSurfAlign_grid3d(grid3d_t *g, std::vector<LevelsetCode> &ls);

    virtual ~MultiSurfAlign_grid3d()
    {
    }

    void clear_patches();
    void add_patch(unsigned code, const char *name);
    void compile_patches();
    std::string patches_to_string() const;

    // Inherited from MeshDomain; overrides implementation from SurfAlign_grid3d
    virtual void closest_point(Mesh::EntityHandle handle,
                               const Vector3D& position, Vector3D& closest,
                               Vector3D& normal, MsqError& err) const;

    virtual void domain_DoF(const Mesh::EntityHandle* handle_array,
                            unsigned short* dof_array,
                            size_t num_handles, MsqError& err) const;

    // methods introduced here

    // instructs the domain that given handle should be projected onto nearest 
    // curve.  also domain_DoF returns 1 for this handle
    void snap_point_to_curve(Mesh::EntityHandle handle,
                             const Vector3D &position);

    // instructs the domain_DoF to return 2 for the given handle.  It will
    // return 3 for any handle that isn't known to the domain.
    void snap_point_to_surface(Mesh::EntityHandle handle,
                               const Vector3D &position);

    // retreives stored codes for given handles; for unknown handles, the
    // returned values are 0.  codes are updated in closest_point(), which
    // must be called prior to get_codes();
    virtual void get_codes(const Mesh::EntityHandle *handles,
                   unsigned codes[], size_t count);

    void set_codes(const Mesh::EntityHandle *handles,
                   const unsigned codes[], size_t count);

    inline bool is_handle_on_curve(Mesh::VertexHandle h)
    {
        return curve_of_handle.find(h) != curve_of_handle.end();
    }

    void clear_handles(void);

    friend class SurfacePatch;
    class SurfacePatch;
    std::vector<SurfacePatch> patches; // elements of dim 2
    friend class SurfaceCurve;
    class SurfaceCurve;
    std::vector<SurfaceCurve> curves; // elements of dim 1
    struct SurfaceVertex;
    std::vector<SurfaceVertex> vertices; // elements of dim 0

private:
    // we need to modify these in closest_point, so mutable:
    mutable std::map<Mesh::VertexHandle, int> curve_of_handle;
    mutable std::map<Mesh::VertexHandle, int> patch_of_handle;

protected:
    virtual size_t find_closest_node
    (Mesh::EntityHandle handle, const Vector3D position) const;


    std::vector<size_t> patch_id_of_face;
    std::vector<std::set<size_t> > patch_ids_of_node;
    std::vector<std::set<size_t> > patch_ids_of_edge;
    std::vector<size_t> curve_id_of_edge;


    void compile_vertices();
    void compile_closed_curves();
    void push_closed_curve(size_t patch_id);

    size_t find_curve(const SurfaceCurve &c) const;
    size_t get_edges_of_patch(const SurfacePatch &patch,
                              std::vector<size_t> &edge_ids);

    void project_to_curve(size_t sn_id, const Vector3D& position,
                          Vector3D& closest, Vector3D& normal) const;

    //////////////////////////////////////////////////////////////////////
public:

    struct SurfaceVertex
    {
        size_t sn_id;
        std::set<size_t> curve_ids;
        std::set<size_t> patch_ids;

        inline SurfaceVertex(size_t id) : sn_id(id)
        {
        };

        inline bool is_on_curve(size_t c) const
        {
            return curve_ids.find(c) != curve_ids.end();
        }

        inline bool is_on_patch(size_t p) const
        {
            return patch_ids.find(p) != patch_ids.end();
        }

    };
    //////////////////////////////////////////////////////////////////////
public:

    class SurfaceCurve : public AlignDomain
    {
        MultiSurfAlign_grid3d *surf_domain;
    public:
        bool closed;
        unsigned code;
        size_t patch_ids[2];
        size_t vertex_ids[2];
        std::vector<size_t> sn_ids; // indices of nodes in surf_nodes that belong to curve
        std::vector<size_t> se_ids; // indices of edges in surf_edges that belong to curve

        std::vector<Mesh::VertexHandle> handles;

        SurfaceCurve(MultiSurfAlign_grid3d *s_d) : surf_domain(s_d)
        {
        }

        bool operator==(const SurfaceCurve & rhs) const
        {
            return (patch_ids[0] == rhs.patch_ids[0]) &&
                    (patch_ids[1] == rhs.patch_ids[1]) &&
                    (closed == rhs.closed) &&
                    (vertex_ids[0] == rhs.vertex_ids[0]) &&
                    (vertex_ids[1] == rhs.vertex_ids[1]);
        }

        bool find_handle(Mesh::VertexHandle h) const
        {
            bool found =
                    std::binary_search(handles.begin(), handles.end(), h);
            return found;
        }

        size_t find_sn_id(size_t sn_id) const
        {
            std::vector<size_t>::const_iterator found;
            found = std::find(sn_ids.begin(), sn_ids.end(), sn_id);
            if (found == sn_ids.end())
                return (size_t) - 1;
            else
                return found - sn_ids.begin();
        }

        virtual LevelsetCode levelset_code
        (Mesquite::Mesh::EntityHandle entity_handle,
         const Mesquite::Vector3D&coordinate) const;

        virtual void snap_to(Mesh::EntityHandle entity_handle,
                             Vector3D &coordinate) const;
        virtual void normal_at(Mesh::EntityHandle entity_handle,
                               Vector3D &coordinate) const;
        virtual void normal_at(const Mesh::EntityHandle* handles,
                               Vector3D coordinates[],
                               unsigned count,
                               MsqError& err) const;
        virtual void closest_point(Mesh::EntityHandle handle,
                                   const Vector3D& position,
                                   Vector3D& closest,
                                   Vector3D& normal,
                                   MsqError& err) const;
        virtual void domain_DoF(const Mesh::EntityHandle* handle_array,
                                unsigned short* dof_array,
                                size_t num_handles,
                                MsqError& err) const;
    };
    //////////////////////////////////////////////////////////////////////
public:

    class SurfacePatch : public Mesquite::MeshDomain
    {
        MultiSurfAlign_grid3d *surf_domain;
    public:
        unsigned code;
        char name[16];

        std::vector<size_t> sn_ids; // indices of nodes in surf_nodes that belong to patch
        std::vector<size_t> sf_ids; // indices of faces in surf_faces that belong to patch
        std::vector<size_t> se_ids; // indices of edges in surf_edges that belong to patch

        std::set<size_t> vertex_ids;
        std::vector<size_t> curve_ids; // a curve contains the boundary edges and nodes
        std::set<size_t> adj_patch_ids;

        std::vector<Mesh::VertexHandle> handles;

        SurfacePatch(MultiSurfAlign_grid3d *s_d,
                     unsigned c, const char *n = "\0")
        : surf_domain(s_d), code(c)
        {
            strncpy(name, n, 15);
            name[15] = 0;
        }

        inline bool has_code(unsigned c) const
        {
            return (code & c) != 0;
        }

        inline bool is_name(const char *n) const
        {
            return strcmp(name, n) == 0;
        }

        bool find_handle(Mesh::VertexHandle h)
        {
            bool found =
                    std::binary_search(handles.begin(), handles.end(), h);
            return found;
        }

        inline std::string to_string() const
        {
            std::ostringstream o;
            std::bitset < 8 * sizeof (code) > bits(code);
            o << (strlen(name) > 0 ? name : "patch") << " : ";
            for (unsigned b = 0; b < bits.size(); b++)
                if (bits.test(b))
                    o << std::setw(3) << b << ", ";
            o << "\b\b;";
            return o.str();
        }

        virtual void snap_to(Mesh::EntityHandle entity_handle,
                             Vector3D &coordinate) const;
        virtual void normal_at(Mesh::EntityHandle entity_handle,
                               Vector3D &coordinate) const;
        virtual void normal_at(const Mesh::EntityHandle* handles,
                               Vector3D coordinates[],
                               unsigned count,
                               MsqError& err) const;
        virtual void closest_point(Mesh::EntityHandle handle,
                                   const Vector3D& position,
                                   Vector3D& closest,
                                   Vector3D& normal,
                                   MsqError& err) const;
        virtual void domain_DoF(const Mesh::EntityHandle* handle_array,
                                unsigned short* dof_array,
                                size_t num_handles,
                                MsqError& err) const;
    };
private:
    void check_codes() const;
    void check_patches_overlap() const;

};

inline void MultiSurfAlign_grid3d::clear_patches()
{
    patches.clear();
}

#endif	/* _MULTISURFALIGN_GRID3D_H */

